!************************************************************
! this subrouitne calculate inner point of isentropic
! supersonic flow
!************************************************************
! Created by : B. G.
! Date       : 2015-08-11
! Revised    :
!************************************************************
! Nomenclature:
!   yp~a00 temporary properties, p-c+, m-c-, 0-streanline
!   lp~T02  coefficients of equations
!   *4f     previous loop properties of point 4
!   *3f     previous loop properties of point 3
!   i,j     loop indicators
!   
!************************************************************
! Warning :
!   1. angle in rad, not degree
!************************************************************

subroutine InnerP()
   ! 1
   !        4()
   !  (3)
   !    2
   use VariableDef
   implicit none
   ! given variables
   ! real*8::x1,y1,V1,theta1,p1,rho1,x2,y2,V2,theta2,p2,rho2
   ! temporary variables
   real*8::yp,pp,rhop,Vp,alphap,thetap,ap,Mp,ym,pm,rhom,Vm,alpham,thetam,am,Mm
   real*8::p00,V0,rho0,a00
   real*8::lp,Qp,Sp,Tp,lm,Qm,Sm,Tm,l0,l12,R0,A0,T01,T02
   real*8::x4f,y4f,p4f,rho4f,V4f,theta4f,x3f,y3f
   integer::i,j
   ! calculate varialbes
   !real*8::x4,y4,V4,theta4,p4,rho4,T4,M4
   
   ! slope of line 12
   l12=(y1-y2)/(x1-x2)
   
   ! set initial c+ and c- values
   yp=y2
   pp=p2
   rhop=rho2
   Vp=V2
   thetap=theta2
   
   ym=y1
   pm=p1
   rhom=rho1
   Vm=V1
   thetam=theta1
   
   theta3=(theta1+theta2)/2
   theta4=theta3
   ! set front loop values
   x4f=0.0
   y4f=0.0
   p4f=0.0
   rho4f=0.0
   V4f=0.0
   theta4f=0.0
   
   do i=1,icor
      call Thermo(pp,rhop,Vp,ap,T4,Mp)
      call Thermo(pm,rhom,Vm,am,T4,Mm)
      alphap=asin(1/Mp)
      alpham=asin(1/Mm)
      ! calculate corfficients
      lm=tan(thetam-alpham)
      Qm=sqrt(Mm**2-1)/(rhom*Vm**2)
      Sm=delta*sin(thetam)/(ym*Mm*cos(thetam-alpham))
      lp=tan(thetap+alphap)
      Qp=sqrt(Mp**2-1)/(rhop*Vp**2)
      if (yp==0) then
         Sp=delta*sin(thetam)/(ym*Mp*cos(thetap+alphap))
      else
         Sp=delta*sin(thetap)/(yp*Mp*cos(thetap+alphap))
      end if
      ! calculate x4, y4
      x4=(-lm*x1+lp*x2+y1-y2)/(lp-lm)
      y4=(-lm*lp*x1+lm*lp*x2+lp*y1-lm*y2)/(lp-lm)
      if (y4 < 0) then   ! if y4 < 0, return
         flagy40=1
         return
      end if
      ! calculate T+, T-
      Tp=-Sp*(x4-x2)+Qp*p2+theta2
      Tm=-Sm*(x4-x1)+Qm*p1-theta1
      ! location of point 3
      x3f=0.0
      y3f=0.0
      do j=1,icor
         theta3=(theta3+theta4)/2
         l0=tan(theta3)
         ! determine x3, y3, theta3
         x3=(-l12*x2+l0*x4+y2-y4)/(l0-l12)
         y3=(-l0*l12*x2+l0*l12*x4+l0*y2-l12*y4)/(l0-l12)
         theta3=theta2+(theta1-theta2)*(y3-y2)/(y1-y2)
         ! justify and exit
         if ((abs(x3-x3f) <= 1e-6) .and. (abs(y3-y3f) <= 1e-6)) exit
         if (i==1) theta4=theta3
         x3f=x3
         y3f=y3
      end do
      ! interpolate for V3, p3 and rho3
      V3=V2+(V1-V2)*(y3-y2)/(y1-y2)
      p3=p2+(p1-p2)*(y3-y2)/(y1-y2)
      rho3=rho2+(rho1-rho2)*(y3-y2)/(y1-y2)
      if (i==1) then
         p00=p3
         V0=V3
         rho0=rho3
      else
         p00=(p3+p4)/2
         V0=(V3+V4)/2
         rho0=(rho3+rho4)/2
      end if
      call Thermo(p00,rho0,V0,a00,T4,M4)
      ! determine R0, A0, T01, T02
      R0=rho0*V0
      A0=a00**2
      T01=R0*V3+p3
      T02=p3-A0*rho3
      ! determine p4, theta4, V4, rho4
      p4=(Tm+Tp)/(Qp+Qm)
      theta4=(-Qp*Tm+Qm*Tp)/(Qp+Qm)
      V4=(-p4+T01)/R0
      rho4=(p4-T02)/A0
      ! justify if point 4 arrived
      if ((abs(x4-x4f) <= el) .and. (abs(y4-y4f) <= el) .and.&
         &(abs(p4-p4f) <= ep*p4f) .and. (abs(rho4-rho4f) <= erho*rho4f) .and.&
         &(abs(V4-V4f) <= eV*V4f) .and. (abs(theta4-theta4f) <= etheta*theta4f)) exit
      
      ! reset c+ and c- values
      yp=(y2+y4)/2
      pp=(p2+p4)/2
      rhop=(rho2+rho4)/2
      Vp=(V2+V4)/2
      thetap=(theta2+theta4)/2
      ym=(y1+y4)/2
      pm=(p1+p4)/2
      rhom=(rho1+rho4)/2
      Vm=(V1+V4)/2
      thetam=(theta1+theta4)/2
      ! store as front loop values
      x4f=x4
      y4f=y4
      p4f=p4
      rho4f=rho4
      V4f=V4
      theta4f=theta4
   end do
   call Thermo(p4,rho4,V4,ap,T4,M4)
   Vv4=V4*sin(theta4)
   Vu4=V4*cos(theta4)   
!   c-line cross detect
   if (x4 <= x2) then
      flagshock=1
   end if
   
end subroutine InnerP
